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Abstract 

An orbiting black hole binary will generate strong gravitational radiation signatures, making 
these binaries important candidates for detection in gravitational wave observatories. The gravita- 

m 

tional radiation is characterized by the orbital parameters, including the frequency and separation 
at the inner-most stable circular orbit (ISCO). One approach to estimating these parameters relies 
on a sequence of initial data slices that attempt to capture the physics of the inspiral. Using cal- 



culations of the binding energy, several authors have estimated the ISCO parameters using initial 
data constructed with various algorithms. In this paper we examine this problem using confor- 
mally Kerr-Schild initial data. We present convergence results for our initial data solutions, and 
give data from numerical solutions of the constraint equations representing a range of physical 

o ■ 

■ configurations. In a first attempt to understand the physical content of the initial data, we find 

o . 

CO , that the Newtonian binding energy is contained in the superposed Kerr-Schild background before 



o 
cr 



the constraints are solved. We examine some deficiencies with the initial data approach to orbiting 
binaries, especially touching on the effects of prior motion and spin-orbital coupling of the angular 
momenta. Making rough estimates of these effects, we find that they are not insignificant com- 
pared to the binding energy, leaving some doubt of the utility of using initial data to predict ISCO 



X 

parameters. In computations of specific initial-data configurations we find spin-specific effects that 
are consistent with analytical estimates. 
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I. INTRODUCTION 



The computation of gravitational wave production from the interaction and merger of 
compact astrophysical objects is a challenge which, when solved, will provide a predictive 
and analytical resource for the upcoming gravitational wave detectors. A binary black hole 
system is expected to be the strongest possible astrophysical gravitational wave source. 
In particular, one expects a binary black hole system to progress through a series of quasi- 
equilibrium states of narrowing circular orbits as it emits gravitational radiation. In the final 
moments of stellar mass black hole inspiral, the radiation will be detectable in the current 
(LIGO-class) detectors. If the total binary mass is of the order of 10M Q , the moment of final 
plunge to coalescence will emit a signal detectable by the current generation of detectors 
from very distant (Gpc) sources. 

Detecting gravitational radiation is also a significant technical challenge. Gravitational 
waves couple very weakly to matter, and the expected signals are much smaller in amplitude 
than ambient environmental and thermal noise. The successful detection of these waves, 
therefore, requires some knowledge of what to look for. In this regard, an orbiting binary 
black hole system is an ideal candidate for detection since the orbital motion produces regular 
gravitational radiation patterns. In such an inspiraling black hole system, the strongest 
waves are emitted during the last several orbits, as the holes reach the innermost quasi- 
stable orbit (here abbreviated ISCO), and as they continue through the final plunge. The 
dynamics of the holes during these final orbits, especially the orbital angular velocity, cuisco, 
and separation, ^iscch determine the dominant characteristics of the detectable waves. Any 
knowledge of these parameters is advantageous for detecting radiation from these binary 
systems. 

The proper way to predict gravitational waveforms for orbiting black holes is to set initial 
data for two widely separated holes, and then solve the evolution equations to follow the 
inspiral through merger and beyond. This problem is well beyond the capabilities of current 
evolution codes. Therefore, to obtain some information about orbiting black holes we, and 
others , turn to the initial value problem. For an introduction to the 

literature, see the review by Cook |9(. Given a collection of initial data for black holes in 
circular orbits with decreasing radius, one tries to identify a sequence of initial data that 
corresponds to instantaneous images of a time-dependent evolution. Circular orbits are 
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chosen because orbits in the early stages of an inspiral are predicted to become circularized 
because of the stronger gravitational radiation near periapse [3] . When a suitable sequence 
of initial data slices has been obtained, they can then be used to determine various orbital 
parameters. For example, the change in binding energy with respect to the orbital radius 
allows one to identify ^isco> an d a similar analysis of the angular momentum gives cuisco The 
difficulty in this approach comes in ensuring that the initial data at one radius correspond to 
the same physical system as the data for another radius. This can be done for some systems 
by using conserved quantities. For example, in the case of neutron stars, constant baryon 
number is an unambiguous indicator of the sameness of the stars. However, in black hole 
physics it is not available; it is unclear how to determine that two black hole initial data 
sets do, in fact, represent the same physical system. 

The initial data approach to studying binary black holes is thus not without problems. 
These difficulties fall in two broad areas. First, there is no unambiguous way to set initial 
data in general relativity. The current algorithms all require some arbitrary mathematical 
choices to find a solution. For instance, the approach we take requires the definition of 
the topology of a background space and of its metric and the momentum of the metric, 
followed by solution of four coupled elliptic differential equations for variables that adjust 
the background fields. But the choice of the background quantities is arbitrary to a large 
extent. The physical meaning of these mathematical choices is not completely clear, but the 
effect is unmistakable. Data constructed with various algorithms can differ substantially, 
even when attempting to describe the same physical system Q]. The data sets can be 
demonstrated in many circumstances to contain the expected Newtonian binding energy, as 
we show below (i.e., the binding energy of order 0(m 2 /£) agrees with the Newtonian result 
at this order). However, the data can differ significantly at 0(m(m/£) 2 ). These differences 
are attributed to differences in wave content of the data which may reflect possible prior 
motion or may simply be spurious. At present it is neither possible to build prior motion 
into the initial data, nor to specify how radiation is added to the the solution, nor to know 
how much there is. It is known that the circular orbits and the ISCO so determined are 
in fact method-dependent. Furthermore, the methods need not even agree that a specific 
dataset represents a circular orbit; their subsequent evolutions may not agree 

A second problem — and the principal physical difficulty with the initial data method for 
studying black hole binaries — is the lack of unambiguous conserved quantities. The best 
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candidate for an invariant quantity is the event horizon area, Au- This area is unchanging 
for isentropic processes due to the proportionality of with the black hole entropy. One 
can argue that since the quasi-circular orbit is quasi-adiabatic, An is nearly invariant over 
some phase of the inspiral. But the inspiral cannot be completely adiabatic because it cannot 
be made arbitrarily slow; the black holes will absorb an unknown amount of gravitational 
radiation while in orbit and will thereby increase in size. Moreover, the event horizon is a 
global construct of the spacetime, and cannot be determined from a single slice of initial 
data. Therefore, one must use the apparent horizon area, A^n, as an ersatz invariant for 



initial data studies 12L Il3lj . When the hole is approximately stationary, these horizon areas 
may be nearly equal jlJ^In dynamic configurations — as should be appropriate for orbiting 
holes — these horizon areas may differ substantially jisl. 11^. 

We will investigate physical content of initial data, focusing on Kerr-Schild spacetimes. 
We examine binding energy to leading order, and find that in our method of constructing 
the superposed Kerr-Schild data, the background fields contain the Newtonian binding en- 
ergy: the subsequent solution of the elliptic equations yields a only small correction. Using 
numerical solutions we present orbital configurations with solved initial data. We give a 
qualitative discussion of physical effects that may confound any attempt to study inspiral 
via a sequence of initial data, and which may affect the determination of the location of the 
ISCO. We give some computational examples consistent with these qualitative predictions. 

II. REVIEW OF INITIAL DATA CONSTRUCTION IN GENERAL RELATIVITY 

In the computational approach we take a Cauchy formulation (3+1) of the ADM type, 



after Arnowitt, Deser, and Misner VJ. \- 111 such a method the 3-metric g^ is the fundamental 
variable. The 3-metric and its momentum are specified at one initial time on a spacelike 
hypersurface. The ADM metric is 

ds 2 = -(a 2 - fop) dt 2 + 20i dt dx i + gij dx i dx j (1) 

where a is the lapse function and f} 1 is the shift 3-vector. Latin indices run 1, 2, 3 and are 
lowered and raised by g^ and its 3-dimensional inverse g l K a and ft* are gauge functions 
that relate the coordinates on each hypersurface to each other. The extrinsic curvature, 
Kij, plays the role of momentum conjugate to the metric, and describes the embedding of a 
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t = constant hypersurface into the 4-geometry. 

The Einstein field equations contain both hyperbolic evolution equations, and elliptic 
constraint equations. The constraint equations for vacuum in the ADM decomposition are: 

R - KijK ij + K 2 = 0, (2) 

Vj (K ij - g ij K) = 0. (3) 

Here R is the 3-dimensional Ricci scalar, and Vj is the 3-dimensional covariant derivative 
compatible with g^j. These constraint equations guarantee a kind of transversality of the 
momentum (Eq. ©)• Initial data must satisfy these constraint equations; one may not 
freely specify all components of and K^. The initial value problem in general relativ- 
ity thus requires one to consistently identify and separate constrained and freely-specifiable 
parts of the initial data. Methods for making this separation, and solving the constraints as 
an elliptic system, include: the conformal transverse-traceless decomposition Q|; the phys- 



ical transverse-traceless decomposition 



and the conformal thin sandwich decomposition 



which assumes a helical killing vector j^, |20|, |2l| . These methods all involve arbitrary choices 



and do not produce equivalent data. Our solution method uses the conformal transverse- 
traceless decomposition jisj |. 

Solutions of the initial value problem have been addressed in the past by several groups [l|, 
0, 0, 0, ^| . It is the case that until recently, most data have been constructed assuming that 
the initial 3-space is conformally flat. The method most commonly used is the approach 
of Bowen and York j^J, which chooses maximal spatial hypersurfaces and takes the spatial 
3-metric to be conformally flat. This method has been used to find candidate quasi-circular 
orbits by Cook Q|, Baumgarte j^j], and most recently, Pfeiffer et al. j^. 

The chief advantage of the maximal spatial hypersurface approach is numerical simplicity, 
as the choice K = decouples the Hamiltonian constraint from the momentum constraint 
equations. If, besides K — 0, the conformal background is flat Euclidean 3-space, there 
are known that analytically solve the momentum constraint j^. The constraints then 
reduce to one elliptic equation for the conformal factor <p. However, it has been pointed out 
by Garat and Price j3] that there are no conformally flat K = slices of the Kerr spacetime. 
Since we expect astrophysical sources to be rotating, the choice of a conformally flat K = 
background will yield data that necessarily contains some quantity of "junk" gravitational 



radiation. Jansen et al. |2J] have recently shown by comparison with known solutions that 
conformally flat data do indeed contain a significant amount of unphysical gravitational field. 
Another conformally flat K = method recently used by Gourgoulhon, Grandclement, and 
Bonazzola J4 , Q| is a thin sandwich approximation based on the approach of Wilson and 
Mathews (21 1 which assumes the presence of an instantaneous rotation Killing vector to 
define the initial extrinsic curvature. They impose a specific gauge defined by demanding 
that K and the conformal factor remain constant in the rotating frame. Since <j) and K are 
a conjugate pair in the ADM approach, this method solves the four initial value equations 
and one second-order evolution equation. The assumption of a Killing vector suppresses 
radiation or, perhaps more accurately, imposes a condition of equal ingoing and outgoing 
radiation. 

to outline some of the difficulties in finding 



In this paper we use Kerr-Schild data 



the ISCO using the initial data technique. We discuss the extent to which initial data set 
by means of superposed Kerr-Schild black holes limits the extraneous radiation in the data 
sets, and we estimate the accuracy of the extant published ISCO determinations. Recent 
works by Pfeiffer, Cook, and Teukolsky also investigate binary black hole systems using 
Kerr-Schild initial data 8[. 



III. INITIAL DATA via SUPERPOSED KERR-SCHILD BLACK HOLES 



The superposed Kerr-Schild method for setting black hole initial data, developed by 



Matzner, Huq, and Shoemaker [25j, produces data for black holes of arbitrary masses, boosts, 
and spins without relying on any underlying symmetries of any particular configuration. The 
method proceeds in two parts. First, a background metric and background extrinsic curva- 
ture are constructed by superposing individual Kerr-Schild black hole solutions. Then the 
physical data are generated by solving the four coupled constraint equations for corrections 
to the background. Intuitively, the background solution should be very close to the genuine 
solution when the black holes are widely separated, and only small adjustments to the grav- 
itational fields are required to solve the constraints. We show that this is true for large and 
also for small separations. This section briefly reviews the superposed Kerr-Schild method 
for initial data, then gives some analytic results to justify this contention. 
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A. Kerr-Schild data for isolated black holes 

The Kerr-Schild j3] form of a black hole solution describes the spacetime of a single black 
hole with mass, m, and specific angular momentum, a = j/m, in a coordinate system that 
is well behaved at the horizon. (We use uppercase M for calculated masses, e.g., the ADM 
mass, and lowercase m for mass parameters, or when the distinction is not important.) The 
Kerr-Schild metric is 

ds 2 = rj^ dx M dx u + 2H(x a %l„ dx^ dx", (4) 

where rj^ is the metric of flat space, if is a scalar function of x^, and / M is an (ingoing) null 
vector, null with respect to both the background and the full metric, 

rriju = 9^%lu = 0. (5) 

This last condition gives 1°Iq = —l l U. 

The general non-moving black hole metric in Kerr-Schild form (written in Kerr's original 
rectangular coordinates) has 

H = ~2~, 2 2~fl> W 

r z + cr cos^ 

and 

_ / rx + ay ry - ax z\ 

M \ '9i9'9i9')' \ ) 

\ r z + cr r z + a z r J 
where r, 9 (and <fr) are auxiliary spheroidal coordinates, z = r(x, y, z) cos8, and </> is the 
axial angle. r(x,y,z) is obtained from the relation, 

4t4+4=l («) 

r z + a z r A 

giving 



r 2 = l -{p 2 -a 2 ) + J\{p^-a^ + a^\ (9) 



with 



p = yj x 2 + y 2 + z 2_ ( 1Q ) 

Comparing the Kerr-Schild metric with the ADM decomposition Eq. (JIJ, we find that 
the t = constant 3-space metric is: 

g ij = 5 ij + 2Hl i l j , (11) 
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Further, the ADM gauge variables are 



Pi = 2Hl li 



(12) 



and 



a 



1 + 2HII 



(13) 



The extrinsic curvature can be computed from the metric using the ADM evolution 
equation 



K 



1 



13 ~ 2a ^ jl3i + Vi/?j ~~ ^ 



(14) 



where a dot ( ' ) denotes a the partial derivative with respect to time. Each term on the 
right hand side of this equation is known analytically. 



B. Boosted Kerr-Schild black holes 



The Kerr-Schild metric is form-invariant under a boost, making it an ideal metric to 
describe moving black holes. A constant Lorentz transformation (the boost velocity, v, is 
specified with respect to the background Minkowski spacetime) A. a p leaves the 4-metric in 
Kerr-Schild form, with H and transformed in the usual manner: 



x" 3 
H'(x' a ) 

GOO 



H ((A- 1 )^ x'P) , 
K\ I, ((A- 1 )^ x'P) 



(15) 
(16) 
(17) 



Note that 1' is no longer unity. As the initial solution is stationary, the only time dependence 
comes in the motion of the center, and the full metric is stationary with a Killing vector 
reflecting the boost velocity. The solution, therefore, contains no junk radiation, as no 
radiation escapes to infinity during a subsequent evolution. Thus, Kerr-Schild data exactly 
represent a spinning and/or moving single black hole. This is not possible in some other 
approaches, e.g., the conformally flat approach ^24]. 



S 



C. Background data for multiple black holes 

The structure of the Kerr-Schild metric suggests a natural extension for multiple black 
hole spacetimes using the straightforward superposition of flat space and black hole functions 

9ij ~ Vij + 2 iH ik ilj + 2 2 H 2 k ilj + • • • , (18) 

where the preceding subscript numbers the black holes. Note that a simple superposition 
is typically not a genuine solution of the Einstein equations, as it does not satisfy the 
constraints, but it should be "close" to the real solution when the holes are widely separated. 

To generate the background data, we first choose mass and angular momentum parame- 
ters for each hole, and compute the respective H and l a in the appropriate rest frame. These 
quantities are then boosted in the desired direction and offset to the chosen position in the 
computational frame. The computational grid is the center of momentum frame for the two 
holes, making the velocity of the second hole a function of the two masses and the velocity 
of the first hole. Finally, we compute the individual metrics and extrinsic curvatures in the 
coordinate system of the computational domain: 

a9%3 — Vij + 2 aH aU aIj-, (19) 
A Ki m = A g mj (V, A /3i + V, - A9ij) ■ (20) 

Again, the index A labels the black holes. Data for N holes are then constructed in super- 
position 



N 



9ij — Vij + X]2 aB aHaU Alj, 

A 

N 

K = Y.aB A Ki\ 



N 



A 



9n(i 



- -$j) n AKi 



AJ^j) 



(21) 
(22) 
(23) 



A tilde ( ) indicates a background field tensor. The simple superposition of the metric 
from Eq. (|18|) (part of the original specification 2^]) has been modified here with the intro- 



duction of attenuation functions, aB 



28, 



The extrinsic curvature is separated into its 



trace, K, and trace- free parts, Ay, and the indices of Ay are explicitly symmeterized. 
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The attenuation functions represent the physical idea that in the immediate vicinity of 
one hole, the effect of a second hole becomes negligible. Near a black hole the conformal 
background superposition (~) metrics approach the analytic values for the single black hole. 
The attenuation function 2 B (i-B) eliminates the influence of the second (first) black hole 
in the vicinity of the first (second). \B equals unity everywhere except in the vicinity of 
the second black hole, and its first and second derivatives are zero at the singularity of the 
second hole. 

The attenuation function used is 

iS = 1 -exp(-^/2<r 2 ), 
where l\ is the coordinate distance from the center of hole 2, 



l\ = i(p 2 -a 2 ) + ^(p2_ a2 )2 + a V , (25) 

p = - 2 x) 2 + (y- 2y) 2 + (z- izf . (26) 

and a is a parameter. In all examples given in this paper, the masses are equal and a = m 2 . 
Fig. ^ shows a typical attenuation function used in calculating our initial data sets. 

A small volume containing the singularity is masked from the computational domain. 
This volume is specified by choosing a threshold value for the Ricci scalar, typically for 
\R\ > 2/m 2 . For a single Schwarzschild black hole, this gives a spherical mask with a radius 
r ~ 0.73 m. In all cases the masked region lies well within apparent horizons in the solved 
data. In practice we find that a small attenuation region (also inside the apparent horizon) 
is necessary to achieve a smooth solution of the elliptic initial data equations near the mask; 
see Section II. D below. Figures El and El show the Hamiltonian and momentum constraints 
for the background space with and without attenuation. We have not varied the masking 
condition to determine what effect the size of the mask has on the global solution. As 
—ed be.ow, Pfeiffer „ al. have instigated this point Q. 

D. Generating the physical spacetime 

The superposition of Kerr-Schild data described in the previous section does not satisfy 
the constraints, Eqs. and hence are not physical. A physical spacetime can be con- 
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FIG. 1: The attenuation function, \B = 1 — exp(— if /2a 2 ), used to calculate our initial data 
solutions. To indicate the effect of the attenuation function in a binary black hole system, we 
also plot the the background metric function g yy in the vicinity of one hole with and without 
attenuation. The Schwarzschild black holes are placed along the y-axis at ±4m. Here i\ is the 
coordinate distance from the center of the second black hole, and the attenuation function width 



structed by modifying the background fields with new functions such that the constraints 
are satisfied. We adopt the conformal transverse-traceless method of York and collabora- 
to,, Q which consists of a confonnal deposition and a vector potentia, that adjusts 
the longitudinal components of the extrinsic curvature. The constraint equations are then 
solved for these new quantities such that the complete solution fully satisfies the constraints. 

The physical metric, <7y, and the trace- free part of the extrinsic curvature, A^, are related 
to the background fields through a conformal factor 



is a = m 



2 




9ij = <P 9ij, 



(27) 
(28) 
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FIG. 2: The Hamiltonian constraint (units m~ 2 ) calculated for the background space for two 
identical Schwarzschild black holes. The black holes are located on the y-axis at y = ±4 m, and 
have zero initial velocity. The solid curve is the background behavior of the constraint without 
using attenuation functions, and the dashed curve is the constraint with attenuation and a = m 2 . 
The masked region is within the radius r < 0.73m. It can be seen that attenuation does not 
necessarily reduce the constraint, but does smooth it. 



where <fi is the conformal factor, and (Iw) 3 will be used to cancel any possible longitudinal 
contribution to the superposed background extrinsic curvature. w l is a vector potential, and 



(Iwf = VV + VV - -g ij V k w k . 

3 



The trace K is taken to be a given function 



K = K. 



(29) 



(30) 



Writing the Hamiltonian and momentum constraint equations in terms of the quantities in 
Eqs. fl2T jl - (p3n|) . we obtain four coupled elliptic equations for the fields and w % [lsl ]: 



V 2 = (1/8) W +\K 2 <t> 5 - 



- 7 (A ij 
-g ij <fVjK - VjA lj . 



r\A lJ + (lwr)(A lJ + (lw) lJ )), 
2. 
3* 



(31) 
(32) 
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FIG. 3: The y-component of the momentum constraint (units m~ 2 ) calculated for the background 
space of two identical Schwarzschild black holes. The black holes are located on the y-axis at y = 
±4 m, and have zero initial velocity. The solid curve is the background behavior of the constraint 
without using attenuation functions, and the dashed curve is the constraint with attenuation and 
a = m 2 . 

E. Boundary Conditions 



A solution of the elliptic constraint equations requires that boundary data be specified 
on both the outer boundary and the surfaces of the masked regions. This contrasts with 
the hyperbolic evolution equations for which excision can in principle be carried out without 
setting inner boundary data since no information can propagate out of the holes. Boundaries 
in an elliptic system, on the other hand, have an immediate influence on the entire solution 
domain. Using the attenuation functions, we can choose simple conditions, = 1 and 
w l = 0, on the masked regions surrounding the singularities. In practice this inner boundary 
condition is not completely satisfactory because it generates small discontinuities in the 
solution at this boundary. These discontinuities are small relative to the scales in the 
problem, and are contained within the horizon. We have made no attempt to determine 



their global effect on the solution. Pfeiffer et al. 



report a similar observation, and note 
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that the location of the boundary does affect some aspects of the solution, though it has 
little effect on the fractional binding energy or the location of the ISCO. 

The outer boundary conditions are more interesting. Several physical quantities of in- 
terest, e.g., the ADM mass and momenta, are global properties of the spacetime, and are 
calculated on surfaces near the outer boundary of the computational grid. Hence the outer 
boundary conditions must be chosen carefully to obtain the proper physics. We base our 
outer boundary conditions on an asymptotic expansion of the Kerr-Schild metric, which 
relies on the ADM mass and momentum formulas to identify the physically relevant terms 
at the boundaries. We first review these expansions and formulas. 

An asymptotic expansion of the Kerr-Schild metric (p ^> m) gives 

r = p (l + 0(p- 2 )) , (33) 
H = m/p (l + 0(p~ 2 )) , (34) 

k = n l + ^- + 0(p- 2 ), (35) 

where = n % = x % j p. (This is the only place where we do not use the 3-metric to raise 
and lower indices, and n^n 1 = 1). a c is the Kerr spin parameter with a general direction: 
a c = aa c . The shift (Eq. [T2|) is asymptotically 

A = ^(^ + a c ^)+0(p- 3 ). (36) 

The asymptotic expansion of the extrinsic curvature in the stationary Kerr-Schild form (cf. 
Eq. (HU)) is 

2 777, 

aK ab = — (-2n a n b + 5 ab ) 
P 

3m 

-—a {e ajc n b + e bjc n a ) rij 

+^r(n a n b -l5 ab )+0(p- i ). (37) 



p3 y~a-o 3 - 

The terms proportional to a c / p 3 in this expression arise from the transverse components of 
(3 a (f3 a n a = 0); the terms of 0(p -3 ) independent of a c arise from the affine connection. Note 
that a = 1 + 0(p _1 ), and will not affect the ADM estimates below. 

The ADM formulas are evaluated in an asymptotically flat region surrounding the system 
of interest, and in Cartesian coordinates they are 

1 rfd gji 0;i j ^ 



M ^ = wJ W-H dS * (38) 
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PfcADM = h f ( Kki ~ Kbb5k ^ (39) 

^ DM = ^ / {xaK u - x b K ai ) dS\ (40) 

for the mass, linear momentum, and angular momentum of the system respectively I30L [31 1 . 
(All repeated indices are summed.) The mass and linear momentum together constitute a 
4- vector under Lorentz transformations in the asymptotic Minkowski space, and the angular 
momentum depends only on the trace-free components of the extrinsic curvature. 

To compute the ADM mass and momentasfor a single, stationary Kerr-Schild black hole, 
we evaluate the integrals on the surface of a distant sphere. The surface element then 
becomes dS l = n l p 2 dfl, where n % is the outward normal and dQ is the differential solid angle, 
we need to evaluate the metric only to order 0(p _1 ); the differentiation in Eq. (|38|) guarantees 
that terms falling off faster than p~ l do not contribute to the integration. The integrand 
is then ^riip 2 n l dVL and the integration yields the expected ADM mass Madm — m - The 
ADM linear momentum requires only the leading order of of K a b, 0(p~ 2 ); terms falling off 
faster than this do not contribute. The integrand of Eq. (|39|) then becomes — ^n a nbn b p 2 dfl, 
yielding zero for the 3-momentum, as expected for a non-moving black hole. 

At first blush, the integral for the ADM angular momentum Eq. (J4*U|) appears warrant 
some concern: To leading order K a b is 0(p -2 ), and the explicit appearance of x a in the 
integrand suggests that it grows at infinity as O(p), leading to a divergent result. However, 
inserting the leading order term of K a b for a single, stationary Kerr-Schild black hole into the 
integrand of Eq. (|40|). we find that the integrand is identically zero. The 0(p -2 ) terms of K a b 
contain the quantities n a nb and 5 a b, which separately cancel because of the antisymmetric 
form of Eq. (|40|). and a divergent angular momentum is avoided. Including the 0(p~ 3 ) 
terms of K a b, we find J^, DM = e a fe c a c m; the symmetry of the other 0(m 2 p -3 ) terms again 
means they do not contribute. This result for J^ DM thus depends on terms in the integrand 
proportional to a that arise from corresponding terms in j3 l proportional to q l where q a is 
a unit vector transverse to the radial direction, q a n a = 0. Only these terms contribute to 
the angular momentum integral; in particular those terms in (3 l proportional to n l /p do not 
contribute. 

The ADM mass and momenta are Lorentz invariant. For a single, boosted black hole, we 
naturally obtain Madm = I'm and Padm = jrnv. The background spacetime for multiple 
black holes is constructed with a superposition principle, and the ADM quantities are linear 
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in deviations about flat space at infinity. Thus the ADM formula?, evaluated at infinity 
in the superposition, do yield the expected superposition. For example, given two widely 
separated black holes boosted in the x-y plane with spins aligned along the z-axis, we have 

M ADM = i7i^ + 272"i, (41) 

pADM = 0? (42) 

jadm _ ^ ^ imiVl fj + l7ni aj + 27 ( 2 m 2 v 2 b + 2 m 2 a) , (43) 



where \b and 2 b are impact parameters [34J, and the tilde (~) superscript indicates that 
these quantities are calculated with the background tensors g a b and K a b. This superposition 
principle for the ADM quantities in the background data is one advantage of conformal Kerr- 
Schild initial data. (Note, in choosing the center of momentum frame for the computation, 

pADM = q 

is a condition for setting the background data.) 
Consider now the ADM integrals for the solved data. The Hamiltonian constraint be- 
comes an equation for the conformal factor, 0. As this equation is a nonlinear generalization 
of Poisson's equation, asymptotic flatness in the full, solved metric requires that 

0^1 + ^ + O(p" 2 ), (44) 

where C is a (finite) constant. This leads to our outer boundary condition for <p, namely 

d p (p(<f>-l)) | p _oo = 0. (45) 

Furthermore, the linearity of the ADM mass integral gives 

M A dm (solved) — ijim + 2 ^ 2 m + C. (46) 

(Here the absence of a tilde (~) indicates that this mass is calculated using the solved g^.) 
At this point we cannot predict even the sign of C, though \C\ is expected to be small for 
widely separated holes. If \C\ —> oo, then the boundary condition Eq. would fail. The 
existence of solutions using this condition, however, provides evidence that this possibility 
does not occur. 

The boundary condition for w % is more subtle. A priori, we expect w % — >• at infinity, 
but a physically correct solution on a finite domain requires that we understand how w l 
approaches this limit at infinity. We construct our boundary conditions on w k by demanding 
that the ADM angular momentum of the full (solved) system be only finitely different from 
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that of the background (superposed) data. That is, given that {g a b, K a b\ and {g a b, K a b} 
have finite differences at infinity, we demand that J a b — J a b also be finite. Using (Eq. (|28}1) 
and (Eq. (|4()|l). we find for the difference in angular momentum 

Jab - Jab = -^j> (x a V( b Wi) - x 6 V( a Wi)) dS l . (47) 

(0 — ► 1 at infinity, and there is no difference at this order between conformal and physical 
versions of w l and g a b at infinity.) 

We have already evaluated an integral of this form, in the discussion of the Kerr angular 
momentum (see Eq. (|40|) and Eq. (|37j)). where we expressed -K" a & in terms of the Kerr-Schild 
shift vector. In that analysis, we noted that falloff of the form 

Wl -^9± ni + 9l qi + o(p-*), (48) 
P P 

with C\ and C<± constant, and qifi 1 = 0, will give a finite contribution to the angular mo- 
mentum. We therefore take as boundary conditions: 

dpipw'm) = (49) 

d p (p 2 w i {5 i] -n i n j ))=U. (50) 

Figures (0])-© display and w l for a simple configuration. In this case the elliptic 
equations were solved on a domain of ±10 m along each axis with resolution Ax = m/8. 
As can be seen in these figures, the functions and w l actually result in little adjustment 
to the background configurations. Also note that the radial component of w l , w l ni, is the 
dominant function. In the graphs plotted here, which give the functions along the y-axis, 
we find of ||w y ||oo ~ 0.03, while Hm^Hoo ~ 3 x 10~ 3 , and ||0 — 1]]^ ps 0.013. Because of 
the symmetry of the configuration, ||w z ||oo is much smaller. Analytically, w z = on the 
y-axis; computationally we find ||w 2 ||oo ~ 5x 10~ 7 . In fact we find in general that the radial 
component of w % is the dominant function in all directions, consistent with our boundary 
conditions, and consistent with the finding that solution of the constraints has small effect 
on the computed angular momentum. Of course the corrections and w % would be expected 
to be larger, for data describing holes closer together. We show below that this data setting 
method leads to generically smaller corrections than found in other methods, thus allowing 
closer control of the physical content of the data. 
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FIG. 4: <j) along the y-axis connecting two nonspinning holes with orbital angular momentum. The 
holes are boosted in the ± x direction with v = 0.196 and are separated by 10 M. Note that 4> is 
very close to unity everywhere. 
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FIG. 5: w x for the same configuration as in Figure 0J 



IV. BINDING ENERGY IN INITIAL DATA 



As a first step towards understanding the physical content of initial data sets, we examine 
in this section the effect of the presence of a second hole on the horizon areas of a first hole 
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FIG. 6: w y for the same configuration as in Figure 0J w z is numerically zero as expected by 
symmetry. 

and on global features such as the ADM integrals and the binding energy of the pair. 
This analysis is carried out for non-spinning holes to first order in the binding energy. A 
comparison to the Newtonian result indicates that the Kerr-Schild background superposition 
data contain the appropriate physical information at this level. We then consider possible 
spin-related phenomena, estimate their magnitude, and discuss their possible effect near the 
ISCO. 

A. Binding energy in Brill-Lindquist data 

Before discussing the conformal Kerr-Schild data, we first consider Brill-Lindquist data 
for two non-moving Schwarzschild black holes These data are conformally flat, and 

K a b = 0. The momentum constraints are trivially satisfied, and the Hamiltonian constraint 
is solved for a conformal factor: = 1+ m/(2r) + m'/(2r'). Here the two mass parameters 
are m, and m', and r and r' are the distances in the flat background from the holes m and 
m' . 

We find that the apparent horizon areas in the solved data correspond to 

777 77?' 

M AH + M AH = m + m' + —- + 0(r 2 ). (51) 
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The subscript "AH" indicates masses computed from apparent horizon areas, and the sep- 
aration in the flat background space is I [36J. We assume that this mass (computed from 
apparent horizons) is close to the total intrinsic mass of the black holes (which is given by a 
knowledge of the spin — here zero — and the area of the event horizon). The binding energy, 
Ef,, can be computed as the difference of the ADM mass observed at infinity and the sum 
of the horizon masses: 

E b = Madm - M AH - M' AH . (52) 
For Brill-Lindquist data Madm = m + m', so that 

£ . = _fflmrf + 0(Oi m 

which is the Newtonian result. 



B. Binding Energy in Superposed Kerr-Schild Data 

We now calculate the binding energy in superposed Kerr-Schild data (set according to 
our conformal transverse-traceless approach) for a non-moving Schwarzschild black hole at 
the origin, and a second such hole at coordinate distance £ away. (£ is measured in the 
flat space associated with the data construction.) We compute the area of the hole at the 
origin to first order and find that the Newtonian binding energy already appears in the the 
background data prior to solving the constraints. Thus, we have an argument justifying the 
result noted at the end of Section IIII El solving the elliptic constraint equations leads to 
small corrections to the Kerr-Schild background data. 

Let both holes be placed on the z-axis; the first hole with mass parameter m at the 
origin, and a second hole with mass parameter m' at z — £. The holes are well separated, 
and we expand all quantities about the origin in powers of e = m'/£ with e< 1. Using 
Schwarzschild coordinates labeled (r,9,(p) (cf. Eq. (jlJ)-Eq. Q for a = 0), the background 
metric tensor is 

g rr = H \- 2e cos 2 9, (54) 

r 

g r g = — 2er sin ^ cos 6 1 , (55) 
g ee = r 2 + 2er 2 sin 2 9, (56) 
9U = r 2 sin 2 9, (57) 
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with all other components zero. The extrinsic curvature of the second hole, 2K ab , is of 0(e 2 ) 
at the origin, and we have simply K ab = \K ab . Similarly, the trace of the extrinsic curvature 
is K = \K. Finally, the non-zero components of A ab are 

~A„ = (l + 2- + 2ecos 2 ^ , (58) 



where 



C2 

A r e = e — sin 6 cos 6, (59) 
r 

A ee = c 2 (l + 2esin 2 #), (60) 
A^ = c 2 sin 2 6 l . (61) 



° 2 ~ 3 V p + 2M (r + 2m) ' 1 ' 



While K a b is not a function of e, and hence contains no information about the second hole, 
perturbative quantities do appear in A ab . This perturbation in A ab arises because we sum the 
mixed-index components of aA^, and because the full background metric, involving terms 
from each hole, is involved in the symmetrization in Eq. (|23|) . 

To calculate the binding energy we first find the apparent horizon area of the local hole. 
For a single Schwarzschild hole, the horizon is spherical and located at pn = 2m; the area 
of the horizon is 167rm 2 . The effect of the second hole is to distort the horizon along the 
z-axis connecting them, and we define a trial apparent horizon surface as / = 0, where 

f = p-2m-J2aiPi(cos8). (63) 
1 

The expansion of / in Legendre polynomials, Pi, expresses the distortion of the local horizon 
away from the zero-order spherical result. This expansion includes a term describing a 
constant "radial" offset in the position of the apparent horizon, a P . This and the other 
terms defining the surface have the expected magnitude, a\ = O(me). We solve for the 
horizon by placing this expression for the surface into the apparent horizon equation 

+ A ab s a s b - = 0, (64) 

where s l is the unit normal to the trial surface 

f.i 



9 ab Uf,i 



(65) 
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The apparent horizon equation is solved to first order, 0(e). One must evaluate the 
equation at the new (perturbed) horizon location. Let F represent the left-hand side of 
the apparent horizon equation (Eq. (|64j) ) , p = 2m is the horizon surface of the single, 
unperturbed hole, and Ph(#) is the new perturbed horizon. We expand F to first order as 



F( PH (6)) = Fo(po) + ^ 



E a ^=°- (66) 



Pa 



Solving (Eq. ()64|)). the only nonzero coefficients in Eq.()63|) are a = mm'/(3f), ci2 = 
—mm'/{2€). Integrating the determinant of the perturbed metric over the horizon surface, 
p = 2m + J2i a>iPi(cos9), we find the area of the apparent horizon to be 

2t 



A h = 16tt U + ^l + 0(m 2 (m'/£) 2 ), (67) 



corresponding to a horizon mass of Mh = m + mm' j (2£) to Newtonian order, ie to order 
0(e) = Or 1 ). 

In this nonmoving case the total ADM mass is just Madm — m + m! . This leads to the 
Newtonian binding energy at this order 



mm! 



E b = -—. (68) 

Because we work only to lowest order in e, Eq. (|6*K|) had to result in an expression of 
0(me), but it did not have to have a coefficient of unity. Both the conformally flat and 
conformally Kerr-Schild data contain the Newtonian binding energy. However, this result 
is obtained in the superposed background Kerr-Schild metric, while the Brill-Lindquist and 
Cadez data give the correct binding energy only after solving the elliptic constraints. This 
is consistent with the small corrections introduced by and w l (0 ~ 1, \w l \ <C 1) in the 



solved Kerr-Schild data (see Section fill EJ) . This fact — that for a superposed Kerr-Schild 
background the solution of the full elliptic problem modifies the data (and the mass/angular 
momentum computations) only slightly — demonstrates how powerful this choice of data can 
be. 

Furthermore, the Newtonian form of the binding energy (e «C 1) means the correct 
classical total energy is found for orbiting situations. If the holes have nonrelativistic motion, 
their individual masses are changed by order 7 ps 1 + 0(v 2 ) = 1 + 0(e). The binding energy, 
which is already 0(e) and is proportional to the product of the masses, is changed only at 
order 0(e 2 ). The ADM mass, on the other hand, measures 7m, and Madm will be increased 
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by mv 2 /2 (an 0(e) increase) for each hole, leading to the correct Newtonian energetics for 
the orbit. 

The apparent horizon is the only structure available to measure the intrinsic mass of a 
black hole. Complicating this issue is the intrinsic spin of the black hole; the relation is 
between horizon and irreducible mass: 

Ah = I6nm 2 rr = 8nm (^m + \J (m 2 — a 2 ) J (69) 

As Eq. ()69|) shows, the irreducible mass is a function of both the mass and the spin, and 
in general we cannot speci fy th e spin of the black holes. For axisymmetric cases Ashtekar's 
isolated horizon paradigm [lj] gives a way to measure the spin locally. We do not pursue 
the point here since we investigate generic and typically non-axisymmetric situations. 



C. Spin effects in Approximating Inspiral With Initial Data Sequences 

We have seen that the initial data contain the binding energy in a multiple black hole 
spacetime. This information can be used to deduce some characteristics of the orbital 
dynamics, particularly the radius of the circular orbit, £, and the orbital frequency, uo. 
Given a sequence of initial data slices with decreasing separations, we determine Eb for each 
slice. The circular orbit is found where 

fL = °- ™ 

The separation at the ISCO orbit, ^isccb hes at the boundary between binding energy curves 
which have a minimum, and those that do not. The curve for the ISCO has an inflection 
point: 



0. (71) 



The angular frequency is given by 



d£ 2 

^isco = -qj. (72) 

The attempt to model dynamical inspiral seems secure for large separation (£ > 15m), 
though surprises appear even when the holes are very well separated. For instance, Eq. (jfi7j) 
above shows that compared to the bare parameter values, the mass increase is equal for the 
two holes in a dataset. Thus the smaller hole is proportionately more strongly affected than 
the larger one is. 
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The physically measurable quantity in question is the frequency (at infinity) Eq. (J72J) 
associated with the last orbit prior to the plunge, the ISCO. This may be impossible to 
determine by the initial data set method. 

To begin with, isolated black holes form a 2-parameter set (depending on the mass param- 
eter m, and the angular momentum parameter, j = ma). For isolated black holes without 
charge the parameters {m, j} uniquely specify the hole. They are equal, respectively, to the 
physical mass and angular momentum. Every method of constructing multiple black hole 
data assigns parameter values gm and gj to each constituent #(/io/e) in the data set. 

There is substantial ambiguity involving spin and mass in setting the black hole data. 
One must consider the evolutionary development of the black hole area and spin. This is a 
real physical phenomenon which contradicts at some level the usual assumption of invariant 
mass and spin. A related concern arises because it is only the total ADM angular momentum 
that is accessible in the data, whereas one connects to particle motion via the orbital angular 
momentum. 

Consider the behavior of the individual black hole spin and mass in an inspiral. For widely 
separated holes, because the spin effects fall off faster with distance than the dominant mass 
effects do, we expect the spin to be approximately conserved in an inspiral. Therefore it 
should also be constant across the initial data sets representing a given sequence of or- 
bits. But when the holes approach closely, the correct choice of spin parameter becomes 
problematic also. 

Newtonian arguments demonstrate some of the possible spin effects. In every case they 
are a priori small until the orbits approach very closely. However, at estimates for the ISCO, 
the effects begin to be large, and result in ambiguities in setting the data (see Price and 
Whelan 3s|). We will consider these effects in decreasing order of their magnitude. 

For two holes, each of mass m in Newtonian orbit with a total separation of £, the orbital 
frequency is 



From recent work by Pfeiffer et al. |2j, the estimated ISCO frequency is of order mVt = 
0.085, corresponding to £ ~ 6.5 m in this Newtonian approach. 

To compare this frequency, Eq. (|73|). to an intrinsic frequency in the problem, we take the 
lowest (quadrupole) quasi-normal mode of the final merged black hole (of mass 2m) which has 




(73) 
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frequency Imio^ pa 0.37; the quadrupole distortion is excited at twice the orbital frequency. 
(We are using the values for a Schwarzschild black hole in this qualitative analysis.) The 
driving frequency equals the quasi-normal mode frequency when i pa 4 m, as might be 
expected. 

To consider effects linked to the orbital motion on the initial configurations, we can 
first treat the effect of imposing corotation. While we show below that corotation is not 
physically enforced except for very close orbits, it is a fact that certain formulations, for 
instance versions of the "thin sandwich" with a helical Killing vector, require corotation in 
their treatment. For any particular initial orbit, corotation is certainly a possible situation. 

In corotation, then, with Eq. (|73|) . for each hole: 

J = ma — Iuj — 4m 2 (mw). (74) 
The result for the moment of inertia / = 4m 3 is the Schwarzschild value 
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Thus 



a = 4mv / 2(m/£) (3/2) . (75) 
Assume a/m <C 1 , and compute the area of this black hole (27 1 : 



A = 8nm(m + \/m 2 - a 2 ) pa lQirm 2 (l - (a/m) 2 /A). (76) 
The horizon mass computed from this area is 

^/A/(lQir) ^ m(l-A(m/£) 3 ). (77) 

At our estimate of the ISCO orbit, ^isco ~ 6m, this effect is of order of 10% of the 
Newtonian binding energy, distinctly enough to affect the location of the ISCO. (At ^isco ~ 
Qm, a/m pa 0.3 for corotation). 

Two more physical effects are not typically considered in setting data. They are frame 
dragging, and tidal torquing. Within our Newtonian approximations, we will find that these 
effects are small, but not zero as the orbits approach the ISCO. In full nonlinear gravity 
these effects could be substantial precisely at the estimated ISCO. 

The frame dragging is the largest dynamical effect. The orbiting binary possesses a 
net angular momentum. For a rotating mass (here the complete binary system) the frame 
dragging angular rate is estimated as the rotation rate times the gravitational potential at 
the measurement point (27 1. Hence 

mn diag = mu(—)pa( T ) pa— (78) 
25 



This is a/m of order 1% at I = 10m; of order 4% at i = 6m. 

The tidal torquing and dissipative heating of the black holes can be similarly estimated. 
As the two holes spiral together, the tidal distortion from each hole on the other will have 
a frequency which is below, but approaching the quasi-normal frequency. Just as for tidal 
effects in the solar system, there will be lag in the phase angle of the distortion, which we can 
determine because the lowest quasi-normal mode is a dissipative oscillator, driven through 
the tidal effects at twice the orbital frequency: 

q + 2 1 q + ulq = F(u). (79) 

Here m 2 q is the quadrupole moment of the distorted black hole. The parameter 7 is (for a 
Schwarzschild hole of mass 2m)) about 2m-7 = 0.089. In Eq. (J79|) the driving acceleration 
F(oS) is identified with the tidal distortion acceleration. We evaluate it at zero frequency: 

q = F{u = 0)/u 2 (80) 

- p- (81) 
The lagging phase, for driving frequency 2u <C u> , is easily computed to be 

« 4 7 ^ (82) 

= 4?r){z) (83) 

This lagging tidal distortion will produce a tidal torque on the black hole, which we can 
approximate using a combination of Newtonian and black hole ideas. The most substantial 
approximation is that the torque arises from a redistribution of the mass in the "target" 
black hole, of amount Am = mm 2 q = m(m/£) 3 . This mass has separation w 4m. Thus the 
torque on the hole is 

r = sin x (lever arm) x AF (84) 
= sin0 x (4m) x (Am2m 2 /£ 3 ) 
= 8 sin 0m(m/£) 6 
w 32{-f/u ){u/uj )m{m/£) 6 
w 60m(m/£) 15/2 . 

What is most important is the effect of this torque on the angular momentum of the hole 
over the period of time it takes the orbit to shrink from a very large radius. To accomplish 
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this, we use the inspira 



from the orbit; see 
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rate (calculated under the assumption of weak gravitational radiation 
I: 

di 128 fm\ 3 



- 



dt 5 \£J 

Thus 



(85) 



dJ dt , , 

n = T Ti < 86) 

5r /m x 3 



128 
-2m 



(87) 

m\ 9 / 2 , . 

l) • < 88 > 



and 

J(£) « m 2 f ^V /2 ; (89) 



i ) 

assuming that there is minimal mass increase from the associated heating (which we discuss 
just below), this identifies the induced spin parameter a = m(m/£) 7 ^ 2 for an inspiral from 
infinity. 

The estimate a = m(m/£) 7 ^ 2 for an inspiral from infinity assumes the mass of the hole 
has not changed significantly in the inspiral. By considering the detailed behavior of the 
shear induced in the horizon by the tidal perturbation, the growth in the black hole mass 
can be estimated |39|| as 



dm dJ . . 

— = uj—, 90 
dt dt 



leading to a behavior 



m \ 5 



Am(£) w 5m — ; (91) 



Consequently, the change in mass can be ignored until the holes are quite close. However, 
the point is that these Newtonian estimates lead to possible strong effects just where they 
become unreliable, and just where they would affect the ISCO. 



These results are consistent with similar ones of Price and Whelan 



who estimated 



tidal torquing using a derivation due to Teukolsky (^(J. That derivation assumes the 
quadrupole moment in the holes arises from their Kerr character, which predicts specific 
values for the quadrupole moment, as a function of angular momentum parameter a. 

Finally we consider an effect on binding energy shown by Wald and also by Dain. Wald 
directly computes the force for stationary sources with arbitrarily oriented spins. He consid- 
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ered a small black hole as a perturbation in the field of a large hole. The result found 41 1 

was 



mm 



S ■ S' - 3(S • n){S' • h) 
I 3 



(92) 



Here, S, S' are the spin vectors of the sources and n is the unit vector connecting the two 



sources. Dain 42j, using a definition of intrinsic mass that differs from ours, finds binding 
energy which agrees with Wald's Eq. (|9*2*|) at 0(£~ 3 ). This is discussed further in Section 

EH 



V. NUMERICAL RESULTS 



We now turn to computational solutions of the constraint equations to generate physical 
data using the superposed Kerr-Schild data. We first discuss the computational code and 
tests, as well as some of the limitations of the code. Finally, we consider physical conclusions 
that can be drawn from the results. 



Code Performance 



The constraint equations are solved (Eq. ()32j0 with an accelerated SOR solver 43j . The 
solution is iterated until the L 2 norms of the residuals of the fields are less than 10 -10 , far 
below truncation error. Discrete derivatives are approximated with second order, centered 
derivatives. We are limited to fairly small domains, e.g., x % G [—12m, 12m] for a typical m/8 
resolution using 193 3 points. 

To verify the solution of the discrete equations, we have examined the code's convergence 
in some detail. The constraints have known analytical solutions — they should be zero — 
which allows us to determine the code's convergence using a solution calculated at two 
different resolutions. Let S\ be a solution calculated with resolution hi, and S 2 be a solution 
calculated with h 2 , then the convergence factor c\ 2 is 

C=*P (93) 
'°g (fe) 

We constructed a conformal background spacetime with two m = 1 non-spinning black 
holes separated by 6m on the y-axis. The elliptic equations were then solved on grids 
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Convergence (c a &) 
a = m/6 a = m/8 a = m/10 


b = m/8 
b = m/10 
b = m/12 


1.70 

1.77 1.86 

1.79 1.86 1.85 



TABLE I: Convergence data for the Hamiltonian constraint, C , for a solution with two m = 1, 
non-spinning holes at = (0, ±3m, 0) in the conformal background, and outer boundaries at 
x l = ±6m. The solution was calculated at resolutions m/6, m/8, m/10, and m/12. The L 2 
norms of C° were calculated over the entire volume of the domain using a mask of radius lm 
around each hole, while the computational mask has a radius of approximately 0.75 m. This larger 
mask was used to compensate for the slight difference of physical location of the mask at different 
resolutions. The norms are as follows: ||C°(m/6)|| 2 = 0.00389054, ||C°(m/8)|| 2 = 0.00238321, 
||C°(m/10)|| 2 = 0.00157387 and | |C°(m/12)| | 2 = 0.00112328. 

with resolutions of m/6, m/8, m/10 and m/12. Tables IV AHV A| show the convergence 
factors as a function of resolution for the Hamiltonian constraint and the x-component of 
the momentum constraint, C x . The convergence for C y is nearly identical to C x , and as the y- 
axis is an axis of symmetry, C z is identical to C x . Figures [7HI3 show the convergence behavior 
of the constraints along coordinate lines. The constraints calculated at lower resolutions are 
rescaled to the highest resolution by the ratio of resolutions squared. We see second order 
for all components with the exception of the points nearest to the inner boundary. 

Solutions of elliptic equations are well-known to be dependent on all boundary data. 
The outer boundary is an artificial boundary, as the the physical spacetime is unbounded. 
Boundary data for this outer boundary are derived from the asymptotic behavior of a sin- 
gle Kerr black hole. On very large domains these conditions should closely approximate 
the expected field behavior, but on small domains these boundary data may only crudely 
approximate the real solution. This error in the boundary data then contaminates the en- 
tire solution, as expected for elliptic solutions. Additional error arises in the calculation 
of the ADM quantities, as spacetime near the outer boundary does not approach asymp- 
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FIG. 7: The Hamiltonian constraint (units m~ 2 ) along y-axis after solving the elliptic equations 
for 4 different levels of resolution. The constraints are rescaled by the ratio of the resolutions 
squared, showing second order convergence. The two non-spinning, instantaneously stationary 
holes of m = 1 are positioned at ±3 on the y-axis. 
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FIG. 8: y-component of momentum constraint (units m~ 2 ) along the y-axis after solving the elliptic 
equations for 4 different levels of resolution, showing second order convergence. The background 
physical situation is the same as in Figure (|7|). The other momentum constraint components 
evaluated on this axis are zero by symmetry, both analytically, and computationally (0(1CP 12 )). 
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Convergence (c a &) 
a = m/6 a = m/8 a = m/10 


b = m/8 
b = m/10 
b = m/12 


1.93 

1.99 2.06 

1.99 2.03 1.99 



TABLE II: Convergence data for the x-component of the momentum constraint, for the same 
configuration as Table IV Al The norms of C x are as follows: | \C x (m/G)\ I2 = 0.00541231, 
\\C x (m/8)\\ 2 = 0.00310937, ||C x (m/10)||2 = 0.00196156 and ||C x (m/12)||2 = 0.00136514. Con- 
vergence factors were also calculated for C y and C z , and found to be essentially identical to the 
data shown here, and thus are not given separately. 
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FIG. 9: z-component of momentum constraint (units m -2 ) along the z-axis after solving the elliptic 
equations for 4 different levels of resolution, showing second order convergence. Other components 
of the momentum constraint evaluated along this line are zero by symmetry, both analytically and 
computationally (O(10 -12 )). The background physical situation is the same as in Figure (|7jl. The 
behavior of the x-momentum constraint along the x-axis is identical to this Figure, as required by 
the symmetry of the problem. 



(m/6) x (6/1 2) 2 

(m/8) x (8/1 2) 2 
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Domain 


Madm 


± 8 m 


1.942 m 


± 10 m 


1.964 m 


± 11 m 


1.974 m 


± 12 m 


1.980 m 



TABLE III: Total ADM Mass for two instantaneously stationary, non-spinning holes separated by 
6m on a grid of discretization Ax = m/8 for four different domain sizes. 

totic flatness. As an indication of the error associated with the artificial outer boundaries, 
we calculated solutions with the same physical parameters on grids of differing sizes. The 
boundary effects in the Madm are given in Table 1111} and Fig. El shows a contour plot of 4> 
for equal mass, nonspinning, instantaneously stationary black holes with the outer bound- 
aries at x % = ±12m. As a further demonstration of boundary effects in our solutions, Fig.lTTI 
shows <p for a configuration examined by Pfeiffer et al. jsj]. Their solution, shown in Fig. 8 
of sj], was computed on a much larger domain via a spectral method j^J. Thus, while we 
achieve reasonable results, it is important to remember that the boundary effects may be 
significant. Moreover, we have only considered the effect of outer boundaries, while errors 
arising from the approximate inner boundary condition have not been examined. 

Other derived quantities also show convergence: Fig. H~2l shows the ADM mass Madm for 
two nonspinning black holes at 6m separation, and different resolutions. The fit is 



M 



ADM 



1.941 + 0.067 



m 



0.422 ( | 

m 



(94) 



showing good second order convergence. The angular momentum calculation is less robust, 
but exhibits approximately first order convergence. The fit to Figure which shows J ADM 
for two nonspinning holes with orbital angular momentum, gives: 



J* DM = ( 1.837 - 0.121 ( -f ) + 0.237 ( ) ) e 12 m 2 . (95) 
Compare this to the angular momentum computed for the background: J^ DM = 2.0m 2 . 



1.837-0.121 (—) +0.237 (—) )e 12 m 2 . 
V m J \ m 
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FIG. 10: Contour plot of 4> for two instantaneously stationary, non-spinning holes of mass pa- 
rameter m = 1. The holes are separated by 6 m along the y-axis. The bold circles indicate the 
apparent horizons. 
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FIG. 11: Contour plot representing the same configuration as Fig.Elbut with the holes separated 
by 10 m along the y-axis. Compare to Figure 8 in ref. ^| 
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FIG. 12: The total ADM energy for two momentarily stationary non-spinning black holes separated 
by 6m at various resolutions. The results exhibit second order convergence. 

B. Physics Results 



The small computational domain does not negate the utility of these solutions as initial 
data for the time-dependent Einstein equations. For instance, Figures El and El show 
data for grazing and elliptical orbits. They are currently being incorporated into the Texas 
binary black hole evolution code. While the small domains do mean that that our data do 
not represent the best asymptotically flat results available from this method, we can still 
verify some of the qualitative analytical predictions of the previous section. In particular, 
Figures El El and El show the conformal factor <fi for holes instantaneously at rest at a 
separation of 6m. In Figure El they are non-spinning; in Figures El and El each has Kerr 
parameter a = 0.5m. In one case (Fig. m]) the spins are aligned; in the other (Fig. 1T51) they 
are antialigned. Table IIVI gives the values of the apparent horizon area of each hole, the 
ADM mass, and the binding energy fraction for these configurations. The binding energy is 



4l| in Section HVCl 



consistent with the analytic estimates of Wald 

Wald's computation of the binding energy for spinning holes, Eq. ()92|). gives for parallel 
or antiparallel spins orthogonal to the separation (as in our computational models): 
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m/6 m/8 m/10 m/12 

Resolution 

FIG. 13: The total ADM angular momentum for two non-spinning holes boosted in the ± x direc- 
tion with v = 0.3162 and separated by 6m at various resolutions (background angular momentum 
J& DM = 2.0m 2 ). 




FIG. 14: Conformal factor (j) for two instantaneously stationary holes separated by 6m with spin 
parameter a = 0.5. The spins are parallel and pointed out of the page. Compare to Fig. 1101 Also 
notice the boundary effect on the outermost contour, labeled 0.999. 
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FIG. 15: Conformal factor <p for the same configuration as Fig. I14l except the spins are anti-parallel: 
the the spin of the hole at (0, —3, 0) points into the page. 




« mis im) 



FIG. 16: (/) f° r a grazing collision between two equal mass, non-spinning holes. The holes are 
centered at y = ±lm and boosted toward each other with v x = =f0.5c, respectively. 



^ / mm' S ■ S' , 

E >-[— + —) < 96 > 

with oppositely directed spins showing less binding energy. For our S = 0.5 m 2 , £ = 6 m 
configuration, this is a change of order 0(2 x 10~ 3 ) between the parallel and the antiparallel 
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FIG. 17: <j) for two non-spinning holes boosted perpendicularly to their separation. The holes 
are separated by 10m and boosted with v x = ±0.196, giving the system a background angular 
momentum of Jf 2 DM = 2.0m 2 . The calculated jf 2 DM = 1.91m 2 and M ADM = 1.970m 2 . The 
Newtonian data correspond to an elliptic orbit at apastron. 





4 f 

AH 


Madm 


M AH 


Binding Energy 


Parallel Spin 


53.20 


1.973 


1.065 


-0.157 = -0.147 xM AH 


Antiparallel Spin 


53.17 


1.974 


1.065 


-0.156 = -0.146 xM AH 


Zero Spin 


57.10 


1.980 


1.066 


-0.151 = -0.142 xM AH 



TABLE IV: MadMi Aah and associated quantities calculated for two holes with m = 1.0 on a 
grid (24 m) 3 with resolution Ax = m/8. 
Quantity for a single hole. 

cases. For the spinning cases we compute a change in binding energy between parallel and 
antiparallel spins of roughly half that, with the correct sign. This rough correspondence 
to the analytic result is suggestive. However, the nonspinning case deviates from the ex- 
pectation that its binding energy should be between that of the spinning cases. Based on 
the scatter in the binding energies shown here, we estimate that we have achieved about 
3% accuracy in the binding energy. With the accuracy of our solution and the size of our 
domain, we are unable to present a clearer dependence of binding energy on spin. 
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VI. OUTLOOK 



To an extent, the difficulty in setting data will become less relevant, as good evolutions 
are eventually achieved. Then data can be set for initial configurations with very large 
separation, and the subsequent evolution will tell us the future dynamics. In the shorter 



45| will give us an indication of the 



term, the iBBH program of Thorne and collaborators 
evolution of the black hole parameters in the inspiral, and will allow a closer identification 
of the corresponding initial data sequence. To point out a couple of additional physical 
effects, note that, besides the historical component associated with a lagging tidal distortion, 
there is the familiar fact that most data setting methods are incapable of accounting for 
the previously emitted gravitational radiation. One can then expect that data describing 
hyperbolic encounters will be more accurate than data sets describing circular motion. This 
is because, in hyperbolic encounters, which are set as distant initial configurations, the 
radiation is more planar, and confined to near each hole. The radiation should then both 
be better defined, and should have less effect on the subsequent evolution than in the more 
distorted orbiting data set. In any case, the understanding of these problems is extremely 
significant in understanding the physical content of the configurations we must solve to 
provide waveforms for the new generation of gravitational wave detectors. For more accurate 
computational results, we are undertaking a both multigrid approach [46]. and a spectral 
approach J, [g] and expect to have extended results soon comparable to those of 8, 44 1. An 
ultimate goal of such a solver is to be able to carry out elliptic solutions at every integration 
time step, to enable fully constrained evolutions. 
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